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Kolmogorov scaling bridges linear hydrodynamic stability and turbulence 
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The way in which kinetic energy is distributed over the multiplicity of inertial (intermediate) scales is a 
fundamental feature of turbulence. According to Kolmogorov's 1941 theory, on the basis of a dimensional anal- 
ysis, the form of the energy spectrum function in this range is the —5/3 spectrum. Experimental evidence has 
accumulated to support this law. Until now, this law has been considered a distinctive part of the nonlinear inter- 
action specific to the turbulence dynamics. We show here that this picture is also present in the linear dynamics 
of three-dimensional stable perturbation waves in the intermediate wavenumber range. Through extensive com- 
putation of the transient life of these waves, in typical shear flows, we can observe that the residual energy they 
have when they leave the transient phase and enter into the final exponential decay shows a spectrum that is very 
close to the —5/3 spectrum. The observation times also show a similar scaling. The scaling depends on the 
wavenumber only, i.e. it is not sensitive to the inclination of the waves to the basic flow, the shape-symmetry of 
the initial condition and the Reynolds number 

PACS numbers: 47.20.-k, 47.35.De, Al.21.Ak 
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In fluid systems, stability and turbulence are the two faces 
of the same coin. The existence of equilibrium: in one case 
laminar, and steady in the mean in the other. The link between 
these two faces is transition 

Unfortunately, or fortunately, depending on the circum- 
stances, turbulence is the rule and not the exception in fluid 
motion. When the energy forcing in the system is sufficiently 
high, transition to turbulence occurs in the short or in the long 
term. In principle, stability and turbulence studies are inti- 
mately connected. In practice, the relevant literature is split 
into two quite distinct fields (however, counterexamples exist, 
see e.g. IlZI, [s]]). The main reason is that the stability can be 
defined physically as the ability of a dynamical system to be 
immune to small disturbances, which necessarily leads to the 
linearization of the mathematical formulation. This is some- 
thing which cannot occur in turbulence where a great number 
of several interacting scales are always active and experimen- 
tally observable. 

However, at any instant, laminar systems host a multiplic- 
ity of scales: the small perturbations which randomly enter 
the system and, in the linear framework, evolve independently 
from each other. Although linearity on one hand allows each 
evolution to be determined singularly, on the other, it should 
be recalled that a large number of perturbations (not even 
bounded, in principle) are present at the same time. In this 
work, we have tried to consider and observe the collective be- 
havior of small perturbations, in particular, those filling the 
intermediate range of wave lengths that the system can host. 
The aim is the understanding and discovering of possible sim- 
ilarities with turbulence behavior 

As an example, in order to understand whether, and to 
what extent, spectral representation can effectively highlight 
the nonlinear interaction that occurs among different scales, it 
could be useful to consider the state that precedes the onset 
of both instability and turbulence in flows. In this condition, 
even if stable, the system is however subject to a swarming 
of small arbitrary three-dimensional perturbations that consti- 



tutes a system of multiple spatial and temporal scales subject 
to all the processes included in the Navier-Stokes equations: 
linearized convective transport, linearized vortical stretching 
and tilting, and molecular diffusion. If we leave nonlinear 
interaction of the different scales, the other features are tanta- 
mount to the features of the turbulent state. 

If it were possible to observe such a system, by computing 
and comparing a large set of three dimensional waves, and 
build spectra, it would be possible, among others, to deter- 
mine if a power scaling in the intermediate range exists and, 
in case, to compare it with the exponent of the corresponding 
developed turbulent state (notoriously equal to - 5/3). In the 
case a power scaling exists, two possible situations can there- 
fore appear A - The exponent difference is large, and as such, 
is a quantitative measure of the nonlinear interaction in spec- 
tral terms. B - The difference is small. This would indicate 
a higher level of universality on the value of the exponent of 
the intermediate range (the inertial range in turbulence), not 
necessarily associated to the nonlinear interaction. 

For this purpose, by solving a large number of initial-value 
problems, we have determined a large set of transient solu- 
tions (a database of 2400 solutions, see the Supplemental Ma- 
terial, section 5) for two typical shear flows: a plane channel 
flow, archetype of wall flows, and a two-dimensional bluff- 
body wake, archetype of free-flows, see Fig. 1 a. Pertur- 
bations that randomly arrive in the system undergo a tran- 
sient evolution which is ruled by the initial-value problem as- 
sociated to the Navier-Stokes linearized formulation Hill. 
These problems must be parameterized through the principal 
physical and geometrical quantities that can influence the life 
of perturbation waves, the angle of obliquity, the symmetry, 
the polar wavenumber and Reynolds number For instance, 
the wavelength of the waves can be varied in a range as large 
as the range of scales that typically fill the field when the sys- 
tem is in the so-called fully developed turbulent state. 

There exists many kinds of transient behavior, very differ- 
ent and not all of which is trivial (for a description, see the 
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text below and in particular the overview presented in figure 
2). The transient lives can last a few basic time scales (ex- 
ternal length referred to the velocity scale of the basic flow) 
as well as order lO** time scales. During these lives energy 
can be accumulated, then smoothly lost or lost and acquired 
again. Usually, different inner time scales appear For in- 
stance, the pulsation can change in a discontinuous way before 
the asymptotic state is reached. This very rich scenario is met 
out of any self-interaction and interaction with other waves 
and, to some degree, is reminiscent, at least qualitatively, of 
the turbulence phenomenology. 

So, the question arises: how to compare spectrally the set 
of very different transients of large, intermediate and short 
waves? Our answer is that whatever be the difference in the 
wave lives, a common phase exits: the time interval where 
the wave exits the transient and enters the exponential asymp- 
totic state. We thus consider the residual kinetic energy owned 
by the waves in this interval and build a spectrum in the 
wavenumber space. 

Let us now describe synthetically the basic flows we have 
here analyzed, as well as the relevant transient computation. 
The basic flows for a channel flow are represented by the 
Poiseuille solution (see Fig. lb and SM, section 4), and for a 
wake by the first two order terms of the Navier-Stokes asymp- 
totic solution described in 11211 (see Fig. Ic and SM, section 
4). The channel flow is homogeneous along the streamwise 
and spanwise directions {x and z). The profiles only vary 
with the coordinate y. In the case of the wake, the profiles 
also slowly evolve with x. As a consequence, the flow is not 
perfectly parallel and we consider two fixed longitudinal sta- 
tions in the region of the wake where entrainment is present 
1 1311 . xo = 10, 50 {xq is the distance from the body normal- 
ized over the body length). Different Reynolds numbers (the 
dimensionless control parameter that gives the order of magni- 
tude of the ratio between convection and molecular diffusion) 
are considered for each example of flow: subcritical (steady 
laminar solution: Re = 30 for the wake. Re = 500 for the 
channel), supercritical (unstable laminar wake Re — 50, 100 
and turbulent channel Re ~ 10000). An initial- value problem 
(IVP) for small arbitrary three-dimensional vorticity perturba- 
tions imposed on the basic shear flows is then considered. The 
viscous perturbation equations are combined in terms of the 
vorticity and velocity ll4ll . and are solved by means of a com- 
bined Fourier-Fourier (channel) and Laplace-Fourier (wake) 
transform in t he p lane normal to the basic flow, see SM sec- 
tion 4 and 12, 15, 16]. The exploration is conducted with re- 
spect to physical quantities, such as the polar wavenumber, 
the angle of obliquity, the symmetry of the perturbation, the 
flow control parameter, and, for the wake, which is not paral- 
lel, the position downstream of the body. For further details 
on the formulation and the numerical methods used to solve 
the initial-value problems see the Supplemental Material (SM, 
section 4). 

To measure the temporal evolution of the energy of 
each perturbation, we define the kinetic energy density. 
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FIG. 1. (a) Base flows and perturbation scheme. The flow profiles 
are qualitatively represented in pink. The perturbation is represented 
by the blue and the yellow waves which propagate in the x and z 
directions, respectively, fc = ai + 7k is the wavenumber, is its 
inclination angle with respect to the basic flow, (b)-(c) Initial con- 
ditions and velocity profiles. The red dotted and solid lines are the 
symmetric and asymmetric parts of the initial perturbations, respec- 
tively. The grey and black lines in the wake panel (c) schematically 
show the velocity variation with the Reynolds number (Re) and the 
position in the wake, xq. Grey lines: supercritical Re (the critical 
value is 47), black lines: subcritical value. Dotted lines: far wake 
field, solid lines: intermediate field. 
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yf ais the computational limits of the domain, u, v and w 
are the transformed velocity components of the perturbation. 
We can also define the amplification factor, G, as the kinetic 
energy density normalized with respect to its initial value, 

G(t;a,7) = e(t;a,7)/e(i = 0;a,7). 

In terms of amplification factors, the early transient evolution 
offers very different scenarios for which we present in the fol- 
lowing a summary of relevant cases. For example, for both 
base flows, we have observed that the orthogonal waves are 
always asymptotically stable. However, the perturbations are 
able to reach very high maxima of the amplification factor G 
(of the order of 10^, see Fig. 2b) before the transients are 
extinguished. Non-orthogonal asymmetric waves present as 
well transients which are not trivial at all (see blue and red 
solid lines in Fig. 2a and 2c). These perturbations are slightly 
amplified in the early stage of their lives, then decrease for 
several hundreds of time scales and in the end they grow with 
the same slope of the correspondent symmetric waves (com- 
pare the asymptotic trends of non-orthogonal solid and dotted 
curves in panels a and c of Fig. 2). During the decreasing 
phase, these transients clearly show an initial oscillatory time 
scale associated to a modulation in amplitude of the average 
value of the pulsation in the early transient, and which is dif- 
ferent from the asymptotic value of the pulsation (see the in- 
sets in panels a-c in Fig. 2) ll5ll . Here the pulsation (angular 
frequency), u, is defined as the time derivative of the wave 
phase at a fixed transversal position (see SM, section 4). Thus, 
the system exhibits two distinct temporal oscillatory patterns, 
the first, of transient nature, and, the second one, of asymp- 
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FIG. 2. Collection of transient lives of the perturbations observed 
through the amplification factor G. Left column: the magnitude of 
the wavenumber is fixed, while the obliquity and symmetry vary (see 
the legend on the bottom). The insets (temporal evolution of the 
pulsation, uj) highlight the typical formation of time scales others 
than the flow external time scale and the wave period. Right column: 
transients variation with changes in the wavenumber magnitude. 



totic nature. 

As a general comment, the most important parameters affect- 
ing these configurations are the angle of obliquity, the sym- 
metry, and the polar wavenumber. While the symmetry of the 
disturbance influences the transient behavior to a great extent 
and leaves the asymptotic fate unaltered, a variation in the 
obliquity and in polar wavenumber can significantly change 
the early trend as well as the final stability configuration. 
The asymptotic behavior for the plane wake, for disturbances 
aligned to the flow, is shown to be in excelle nt ag reement 11711 
with 2D spatio-temporal modal analyses il8ill9tl and with the 
laboratory determined frequency and wave length of the par- 
allel vortex shedding at Re = 50 and 100 ll20ll . See also in the 
SM, the section Asymptotic pulsation (section 3, figure S3) 
where information on the channel flow pulsation measured 
in the laboratory and from the present IVP computations are 
given. 

We now come back to the spectral analysis of common 
phase in the lives of the perturbations, that is the transition 
between the end of the transients and the settlement of the 
asymptotic condition. To compare the residual kinetic en- 
ergy of the waves in correspondence to this transition, we 
assumed that when the asymptotic exponential temporal be- 
havior is reached, the temporal growth rate, r, defined as 
r{t; a, 7) = log{e)/{2t) 12}], must approach a real constant. 
In order to determine the temporal region in which the evolu- 
tion behaves exponentially, it is necessary to monitor the in- 
stants beyond which the condition r const is satisfied. Do- 
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FIG. 3. Spectra of the amplification factor G for the collections of 
linear traveling waves observed at the time, Te, where the perturba- 
tions are out of their transient. The spectra obtained in this study have 
been compared qualitatively with the turbulent field spectra available 
in the literature and obtained from laboratory or numerical exper- 
iments. The legend in the bottom panel specifies the symbols as- 
sociated with either the results of the present study or those of the 
laboratory, and numerical experiments carried out on fully turbulent 
flow fields by other authors. 



ing this, we introduce a necessary condition, that, however, is 
not sufficient to determine the instants where the perturbations 
can be compared. In fact, the transition toward the exponential 
behavior is smooth, very long, and different from case to case 
(in some instances, it can be oscillatory). As a result, it is not 
numerically efficient to obtain the instants where r starts to be 
constant. We have thus associated a second condition to be 
satisfied together with the constancy of the growth rate, which 
directly acts on the energy temporal rate of variation. To this 
aim, we have selected the instants at which the amplification 
factor reaches a given rate of variation, either in growth or in 
decay. This situation is represented by the instant, that we call 
observation time, T^, where dG/dt < e or dG/dt > 1/e, 
with e = 10~". Here, n is an arbitrary positive quantity (for 
instance, a positive integer) that we have fixed equal to 4. It 
is possible to show that the present results - in particular, the 
existence of an intermediate spectral range where the spectral 
decay exponent is very close to that of the Kolmogorov theory 
- do not depend on the choice of n (see SM, section 1). 

In Fig. 3 and 4, the results of the measure of the energy 
residuals out of the transients are shown. The spectral val- 
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ues of G, for both the channel flow case (panels a,c) and the 
wake flow case (panels b,d) show a scaling in the intermediate 
range of the polar wavenumber {k G [10, 200] for the channel, 
k € [2, 50] for the wake) that is amazingly close to the turbu- 
lent canonical value of —5/3. For shorter wavelengths, char- 
acterized by very short transients, the scaling is a little higher 
in magnitude, approximatively equal to —2. This result does 
not appear to be influenced to any great extent by the wave 
obliquity, the symmetry, or the Re. However, it is possible 
to observe that purely orthogonal waves show a closer scal- 
ing to —2 than to —5/3, even at intermediate wavenumbers. 
In general, a full decade of intermediate wavenumbers can be 
observed for both the wall flow and the free flow. These data 
gather all the stable waves occurring in the intermediate range 
and in the dissipative range. We would like to point out that 
the data do not highlight a dependence on Re, the flow control 
parameter For longer waves (fc < 10 and A: < 1 — 2 for the 
plane Poiseuille flow and the bluff-body wake, respectively), 
the results depend on the perturbation inclination, the sym- 
metry of the initial condition, and on the boundary conditions 
(geometry of the system). As expected, they do not reveal any 
universal behavior. In Fig. 3, panels a-b-c-d, experimental 
(laboratory and numerical) measurements 12 1142411 in the tur- 
bulent states have been included for the sake of comparison of 
stable linear perturbations and turbulent scales. 

These results appear strengthened by the fact that even the 
observation times, Tg, present the same scaling, see Fig. 4. 
This outcome is not at all trivial. It is sufficient to consider 
that the observation time includes the transient, and that the 
different kind of transients we observed are very complex and 
can vary in length over 4 orders of magnitude when moving 
across the space of parameters (the wavenumber, the symme- 
try, the obliquity, the flow control parameter. Re, and, in the 
case of a basic evolving spatial flow, the position). 

A possible scenario for fields with a sufficiently high 
Reynolds number to be in a turbulent state, which includes 
the present findings, could be the following. The nonlinear 
interaction works intensively on the long waves that are un- 
stable by blocking the linear growth and by settling their ki- 
netic energy to a level that can be associated with the global 
Reynolds number of the system. The surplus of energy that 
the unstable waves potentially gain during the linear evolu- 
tion is then transferred to shorter waves, that are longitudinal, 
oblique or orthogonal to the basic velocity field. The trans- 
fer is particularly intense for asymptotically stable intermedi- 
ate waves (namely the inertial scales, according to the termi- 
nology used in turbulence). It is possible to suppose that the 
transfer is physically triggered during the phases where the 
amplification factor is maximum, see Figure 2. This process 
continues and the energy in the range of intermediate waves 
reaches values that can be experimentally observed, in the lab- 
oratory, or in the so called Navier-Stokes direct numerical so- 
lutions. Taking into account the present results, it is possible 
to say that the nonlinear interaction distributes the relative en- 
ergy over different wavenumbers in a way that corresponds to 
the relative residual energy each wave has when, after the tran- 
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FIG. 4. Spectrum of the observation time where the wave energy is 
measured and compared. The data refer to the entire collection of 
linear perturbations shown in figure 3. Stable basic flows (top pan- 
els), unstable flows (bottom panels). The plane channel flow is on 
the left, the bluff body wake is on the right. Symmetric (circles) 
and asymmetric (triangles) initial perturbations. Obliquity angles: 
(f} — (black), (fi — n/4: (blue) and (f) = n/2 (cyan). It is possible 
to observe an abundant decade where longitudinal, oblique and or- 
thogonal perturbations present a power-law decay which is close to 
-5/3 (red curves). Larger wavenumbers, which are influenced more 
by diffusion-dissipative processes, steep towards a decay exponent 
of about -2. 



sient, it reaches a common decay threshold. This process is 
less efficient on short waves. For these waves, the maximum 
of the amplification factor becomes progressively smaller as 
the wavenumber magnitude increases (see Figure 2b) and, as 
a consequence, it is more difficult for them to couple with 
longer waves. A rapid fall in energy to levels many orders 
of magnitude lower, and, thus closer to the specific levels of 
their linear evolution, is therefore observed in the experimen- 
tal spectra (see Figure 3). 

In conclusion, we have observed a fingerprint of the Kol- 
mogorov scaling inside the collective behavior of transient in- 
termediate perturbation waves, which always are asymptot- 
ically stable. These new observations are not specific of a 
peculiar kind of flow (wall bounded or free). This can mean 
that this scaling is not only one of the major signatures of the 
turbulence interaction, but it also exists hidden inside the dy- 
namics of linear stable waves, where even the self-interaction 
is absent. Since our observations do not depend on the system 
control parameter (Reynolds number), on the kind of initial 
condition and on geometrical parameters, such as the wave in- 
clination, they could also reveal a new set of structural proper- 
ties of the Navier-Stokes equation solutions. In particular, we 
think that they can be used to build a bridge between the linear 
and the nonlinear interaction in multi-scale systems. Given 
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the high variabiUty of the shape and length of the transient life, 
we consider remarkable that the Kolmogorov scaling does not 
appear only in the energy spectra, but also in the spectra of the 
observation times. 
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1 Wave comparison criterium. Residual energy out of 
the transient. 

Figure SI, sec in particular the caption content, shows that the double criterium used to 
define the observation time Tg, where we measure the energy of the traveling waves, is not 
conditioning the existence and the exponent of the spectral scaling presented in Figs. 3 and 
4 of the main paper, and in Fig. S2 of the following section. 

2 Additional spectral distributions 

In figure S2, the reader can find further data on both the energy and observation time spectra 
in the wake. With regard to Fig. 3, in the main text, one can see how the energy, in particular 
that of the unstable waves, is growing with Re and decreasing with the distance from the 
body. 



1 




10" 



Channel flow 
Re=10000 
Sym ic, k=2, (|)=jr/4 

50^ 



100 



10' 



10-^ 



10" 



10" 



10 



-10 



b |dG/dt|<E (stable waves) 
Q |dG/dt|>1/e (unstable waves) 


O 

^-5/3 


Channel flow 
Re=10000 


O £=10-" 


Sym lo, (|)=ji/4 


E=10"^ 
O E=10~^ 





10 



10 



k 10' 



10 



10^ 



10" 



10" 



^ |dG/dt|>1/8=10'*^ 




10 



-10 



Wake flow, Re=100 
Xq=10, Asym ic, (|)=0, k=0.85 



100 200 ^ 300 400 



10= 

G 

10° 



10 



d 

|dG/dt|<£ (stable waves) 
|dG/dt|>1/e (unstable waves) 

O 6=10-" 
O E=10-= 




10" 



10 k 



10 



10 



Figure S 1; Criteria to define the exit time from the transient part of the perturbation 
life, (a)-(c) channel and wake flows: we consider here G, the amplification factor, i.e the 
energy normalized over the value taken at the initial instant, and r, the temporal growth 
rate r(i;a,7) = log{e)/ {2t), and show their time dependence: - (black curves) |r|, - (blue 
curves) d\r\/dt, - (red curves) \dG/dt\). In order to determine the temporal region in which 
the evolution behaves exponentially, it is necessary to monitor the instants in which the 
condition r — > const is satisfied. Since this transition is smooth, we need the help of a 
further condition to precisely define an instant at which the wave energy can be compared. 
To this aim. we have selected the instants at which the amplification factor reaches a given 
rate of variation cither in growth or in decay. This situation can be represented by the 
instant, that wc call observation time, Tg, where dG/dt < e or dG/dt > 1/e, with e = 10~", 
where n is an arbitrary positive quantity (for instance, a positive integer). Panels (b) and 
(d) in this figure show that the present results do not depend on the choice of n. In this 
example, we have tested the values n = 2,3 and 4. 
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Figure S 2: Additional spectral distributions for the wake flow, (left) Spectra of the ampli- 
fication factors G and (right) spectra of Te. (a)-(d) Re = 30, xo = 10. (b)-(e) Re = 50, 
xo = 10. (c)-(f) Re = 100, Xo = 50. 
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Zebib, 1987, Re=30, 2D DNS data 




Pier, 2002, Re=100, 2D DNS data 




Barkley, 2006, Re=100, 2D DNS data 




Norberg, 1994, Re=100, lab data 




Nishioka & Sato, 1978, Re=100, lab data 


□ 


Williamson, 1989, Re=100, lab data 


□□□□ 


Nishioka et al.,1975, Re=3000, 4000, 5000, 6000, 7000, lab data 




Ito, 1974, Re=3000, 4000, 5000, 6000, 8000, theordata 




Asai & Floryan, 2006, Re=4000, 5000, 6000, lab/theor data 



Figure S 3: Spectrum in the wavenumber space of the pulsation of the waves shown in figure 
3 of the main text. Stable basic Sows (top panels), unstable fiows (bottom panels). The 
plane channel flow is on the left, the bluff body wake is on the right. Symmetric (circles) and 
asymmetric (triangles) initial perturbations. Obliquity angles: (f) = (black) and = 7r/4 
(blue). The legend in the bottom panel specifies the symbols associated with either the 
results of the present study or these of the laboratory, together with numerical experiments 
and theoretical analysis carried out by other authors (high Reynolds number channel flows, 
normally in turbulent configuration, arc kept laminar). 
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3 Asymptotic pulsation 



The spectral distribution of the pulsation is presented for the collection of linear perturbations 
shown in Fig. 3 of the main text. The distribution highlights a proportionality relationship, 
i.e. nondispersive behavior. No experimental results covering the three decades that we can 
observe are available in the literature. However, some data are available for the longer wave 
range, the most unstable, thus the most frequently studied in the field of hydrodynamic 
stability [Tl[2l[3jEl[5l|6l[3|8l[9|, see Fig. S3 (channel: panels a, c; wake: panels b, d). 
The contrast is very good. It should be noted, that only Williamson [1] presented a complete 
datum for the wake where the frequency of the wave is coupled to the value of its wavelength. 
For the other results, the assumption is made that the measured frequency is globally relevant 
for the flow in near critical conditions, although a critical wavenumber is not specified. We can 
observe two abundant decades where longitudinal and oblique waves present almost universal 
behavior, regardless of the symmetry of the initial condition and the Reynolds number: the 
frequency is proportional to the wavenumber (red curves, the exponent is 1). It is possible 
to observe the same behavior for the orthogonal waves (not reported here), but the pulsation 
values are 20 orders of magnitude lower (the values range at around 10~^°). This means 
that orthogonal waves, when in asymptotic conditions - i.e. outside their transients - become 
standing waves. 



4 Methodology 



The continuity and Navier-Stokes equations that describe the perturbed system, subject to 
small three-dimensional disturbances, are linearized and written as 

du ^ dv ^ dw Q 
dx dy dz ' 

du jj'^^ ~9C/ dp 1 
dt dx dy dx Re 

dv ^dv dp 1 

dw ^,dw 9p 1 , 

where {u{x,y, z,t), v{x,y, z,t), w(x^y,z,t)) and p{x,y,z,t) are the perturbation velocity 
components and pressure, respectively. U and dU/dy indicate the base flow profile (under 
the near-parallelism assumption) and its first derivative in the shear direction. For the 
channel flow, the independent spatial variable z is defined from — cio to -foo. the x variable 
from — oo to +oo, and the y from —1 to 1. For the wake flow, z is defined from — oo to +oo, 
x from to -foo, and y from — oo to -|-oo. All the physical quantities are normalized with 
respect to a typical velocity (the free stream velocity, t//, and the centcrline velocity, J7o, for 
the 2D wake and the plane Poiseuille flow, respectively), a characteristic length scale (the 
body diameter, D, and the channel half-width, ft,, for the 2D wake and the plane Poiseuille 
flow, respectively), and the kinematic viscosity, v. 

The base flow of the wake is approximated at an intermediate [xq = 10) and at a far 
longitudinal station (xq 50), through two-dimensional analytical expansion solutions flOj 
of the Navier-Stokcs equations. Assuming that the bluff-body wake is a slowly evolving spatial 
system, the base flow is frozen at each longitudinal station past the body, by using the first 

orders (?? = 0,1) of the expansion solutions [TUj, U{y; xq, Re) = 1 — axQ^^^cxp 
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where a is related to the drag coefficient Co [a = |(i?e/7r)^/^CD(i?e), see [T^), and xo is the 
streamwise longitudinal station. 

The plane channel flow is homogeneous in the x direction and is represented by the 
Poiseuillc solution, U{y) = 1 — y^. 

By combining equations ([T]) to (jH) to eliminate the pressure terms, the perturbed system 
is fully described in terms of velocity and vorticity [llj by the following equations: 



d 


ox J 




dv(PU 


di " 




dx dy^ 


d 


VU— 

dx y 




dv dU 


dt " 


) - 


dz dy 



v^v = r, (5) 
^v^r, (6) 



1 

'We 



V'i„ (7) 



where Uy is the transversal component of the perturbation vorticity. On the basis of kine- 
matics, the following relation holds: 

f = ^ _ ^ (8) 
dx dz 

The physical quantity, F, physically links the perturbation vorticity in the x and z directions 
{uJx and uj^, respectively) to the perturbation velocity field. 

The perturbations are then Fourier transformed in the x and z directions for the channel 
flow. Two real wavenumbers, a and 7, are introduced along the x and z coordinates, respec- 
tively. A combined Laplace-Fourier decomposition is performed for the wake flow in the x 
and z directions. In this case, a complex wavenumber a = Ur + iai can be introduced along 
the X coordiuEite, as well as a real wavenumber, 7, along the z coordinate. The perturbation 
quantities {v^T,ujy) involved in the system dynamics are now indicated as {v^T,Cjy), where 

/ + OO P + OC ^ 
/ fix,y,z,t)e-'''--''''dxdz, (9) 
-oo J —oo 

indicates in the a — 7 phase space the two-dimensional Fourier transform (in the case of the 
channel flow) of a general dependent variable, /, and 

/+00 i<+oo 
/ g(x,y,z,i)e-'"--'^^da:dz, (10) 
00 Jo 

indicates the two-dimensional Laplace-Fourier transform (in the case of the wake flow) of a 
general dependent variable, g. To obtain a finite perturbation kinetic energy, the imaginary 
part, tti, of the Laplace transformed complex longitudinal wavenumber can only assume non- 
negative values and can thus be defined as a spatial damping rate in the streamwise direction. 
In so doing, perturbative waves can spatially decay {ai > 0) or remain constant in amplitude 
{ai = 0). Here, for the sake of simplicity, we have ai = 0, therefore a ~ Ur- The governing 
partial differential equations we consider are thus 



dy'^ 

dr 

'dt 

dCjy 

'df 



k^v = F, 



d'^U 1 ( d'^r A 
-ikcos((j))Ur + ikcos{<j))—-irV -\- — — -— r — fc^r , 
dy^ Re \ dy^ I 



/ ,Nrr- ■ , :xd.U ^ 1 / d-uj. 

-ikcosimU ijjy — iksin[(p)——v H — - 

dy Re \ dy 



(11) 
(12) 

(13) 
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Channel flow 


Wake flow 




or 


n{a,-i)exp{~y'^)cos{y) 
or 

il(a, ^)exp{—y'^)sin{y) 


iOy{y,t = 0) 









Table S 1: Initial conditions for the channel and wake flows. 



where (j) = tan ^{j/a) is the perturbation obliquity angle with respect to the x-y plane, 
k ~ ^ oP- -V is the polar wavenumber and a ~ kcos{(j)), 7 = ksin{(j)) are the wavenumber 
components in the x and z directions, respectively. 

Various initial conditions can be used to explore the transient behavior. The important 
feature here is the ability to make arbitrary specifications. It is physically reasonable to 
assume that the natural issues affecting the initial conditions are the symmetry and the 
spatial lateral distribution of disturbances. It has been observed [T2l[l3l[T4] that, keeping all 
the other parameters fixed, if the perturbation oscillates rapidly or mainly lies outside the 
shear region then, for a stable configuration, the final damping is accelerated while, for an 
unstable configuration, the asymptotic growth is delayed. However, the general qualitative 
scenario is not altered. Therefore, to perform a more synthetic perturbative analysis, we 
only focus on symmetric and asymmetric inputs which are localized and distributed over 
the whole shear region. The transversal vorticity Cjy{y,t) is initially taken equal to zero to 
highlight the three-dimensionality net contribution on its temporal evolution. The effects of 
non-zero initial conditions on the transversal vorticity ujy{y,t) can be found in }13[ I14| . The 
imposed initial conditions are reported in Table SI, for the channel and wake flows. ri(a,7) 
is the phase space transform of the x-z dependence prescribed at time t = 0. Here, we set 
fi(a,7) = 1, which means that no wavenumber is initially biased in the phase space. 
For the channel flow no-slip and impermeability boundary conditions are imposed, 

dv 

v{y = ±l,t) = ^(2/ = ±1>0 = ^yiy = ±1,^) = 0, (14) 
while for the wake flow uniformity at infinity and finiteness of the energy are imposed. 



v[y — ?> ±00, t) = —(y ^ ±00, t) = ujy{y — > ±00, t) 
dy 



0. 



(15) 



To measure the growth of the perturbations, we define the kinetic energy density, e, 

1 r+Vf 



e(t;a,7) 



{]u\^ + \v\^ + \u,\')dy, 



where — ?// and y/ are the computational limits of the domain, while u, v and w are the trans- 
formed velocity components of the perturbation field. We can also define the amplification 
factor, G, as the kinetic energy density normalized with respect to its initial value: 



G(i;a,7) 



e(^;a,7) 
e{t = 0; a, 7) 



(16) 



Assuming that the temporal asymptotic behavior of the linear perturbations is exponen- 
tial, the temporal growth rate, r, can be defined as r{t]a,j) = log{e)/{2t). 
This quantity has a precise meaning when the asymptotic state is reached, that is, when it 
becomes a constant. 
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The angular frequency (pulsation), u, of the perturbation can be defined considering the 
phase, Lp, of the complex wave at a fixed transversal station j/o (for example, j/o = 1 for the 
wake fiow and yo = 0.5 for the channel flow) 




(17) 



(18) 



Although defined at any time t, the frequency lo is here referred to as an asymptotic property 
of the perturbation. 

Equations (|lip - (fT5|) are numerically solved by the method of lines: the equations are first 
discretized in the spatial domain and then integrated in time. The spatial derivatives in the 
y domain are discretized using a second-order finite difference scheme for the first and second 
derivatives. One-sided differences are adopted at the boundaries, while central differenced 
derivatives are used in the remaining part of the domain. The spatial grid is uniform with a 
spatial step, h, which is equal to 0.05 and 0.004, for the wake and channel flows, respectively. 
Since the wake flow is spatially unbounded in the transversal direction, the spatial domain, 
[~2//iZ//]i is chosen so that the numerical solutions are insensitive to further extensions of 
the computational domain size. Equations ()lip -(fT3 |) are then integrated in time by means of 
an adaptative one-step solver, based on an explicit Runge-Kutta 3(2) formula. 

5 Database setting 

The Matlab scripts and the set of solutions of the linearized Navier-Stokes perturbative 
equations can be found at https://130.192.25.166 (username: guest, password: etipso). By 
using the "File Browser" it is possible to access the disks 
/LaCie (usb)#2/Linearized_NS_solutions 

For the channel flow, the perturbative analysis considers 4 parameters (Reynolds number, 
symmetry/ asymmetry, angle of obliquity and wavenumber) and 660 solutions are obtained. 
If one, by using the "File Browser", goes to 
/LaCie (usb)#2/Linearized_NS_solutions/ channel, 

he will find the database organized as in Fig. S4. Every folder, corresponding to a certain 
specification of the above parameters, contains the following text files: 

• prefix_t_n.txt: the temporal points, M, at which the solutions are computed through 
the Matlab code, are reported in column; 

• prefix_u_n.txt, prefix_v_n.txt, prefix_w_n.txt: these files contain the perturba- 
tion velocity field components. Each of them has two columns, the first one for the 
real part and the second one for the imaginary part of the velocity component. The 
column length is MxN, where A'' is the number of spatial grid points (TV = 501 for 
the channel flow) and M are the temporal instants. For every fixed time, the velocity 
spatial distributions are put in column; 

• prefix_omega_y_n.txt: this file contains the transversal vorticity component and is 
structured as the above velocity field files; 

• prefix_energy_n.txt: in this file, the kinetic energy density, e, is reported. The first 
column expresses the temporal points, Af, and the second one the corresponding energy 



values. 
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Figure S 4: Database structure for the channel flow. 
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Figure S 5: Database structure for the wake flow. 
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Each of these files contains, as a prefix in its name, the parameter information and has an 
increasing number, n, as a suffix. This integer number, n, accounts for the fact that outputs 
are periodically saved (after a variable temporal interval) as the equations are integrated in 
time. 

Let us make an example. Suppose we are interested in the simulation with parameters 
Re = 500, symmetric initial conditions, = 7r/4, fc = 3. If one goes to 
/channel/Re_500/Re_500_sym/Re_500_sym_phi_45/Re_500_sym_phi_45_k_3, 
he will find the above text files with prefix Re_500_s3nn_phi_45_k_3 and suffix n > I. 

For the wake flow, data are organized in an analogous way to the one described for the 
channel flow and can be found in the following folder, 
/LaCie(usb)#2/Linearized_NS_solutions/waJs;e. 

It should be recalled that 5 parameters are here considered (Reynolds number, wake position, 
symmetry/asymmetry, angle of obliquity and wavenumber) and 1740 solutions are computed 
(see Fig. S5, where the scheme of the database structure for the wake flow is shown). 
Therefore, in the corresponding folders and files, as part of the suffix, information on the 
chosen downstream station, xq, are added similarly to what done for the other parameters. 
Moreover, concerning the perturbation velocity and the transversal vorticity files, it should 
be mentioned that the spatial domain is enlarged when long perturbative waves are analyzed. 
Therefore, the number of grid points, N, depends on the wavenumber considered: N = 801 
if fc > 1, 1201 if fc G [0.75, 1], iV = 1601 if fc 6 [0.45, 0.7]. 
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